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We present a new analysis of the anisotropic spectral energy distribution in incompressible mag- 
netohydrodynamic (MHD) turbulence permeated by a strong mean magnetic field. The turbulent 
flow is generated by high-resolution pseudo-spectral direct numerical simulations with large-scale 
isotropic forcing. Examining the radial energy distribution for various angles 6 with respect to Bo 
reveals a specific structure which remains hidden when not taking axial symmetry with respect to 
Bq into account. For each direction, starting at the forced large-scales, the spectrum first exhibits an 
amplitude drop around a wavenumber ko which marks the start of a scaling range and goes on up to 
a dissipative wavenumber kd(6). The 3D spectrum for k > ko is described by a single ^-independent 
functional form F(k/kd), the scaling law being the same in every direction. The previous properties 
still hold when increasing the mean field from Bo = 5 up to Bo = 10 b rms , as well as when passing 
from resistive to ideal flows. We conjecture that at fixed Bo the direction-independent scaling regime 
is reached when increasing the Reynolds number above a threshold which raises with increasing Bo. 
Below that threshold critically balanced turbulence is expected. 



It is known that in presence of a mean magnetic field 
assumed here to point in the z-direction, Bo = BqG Zi 
nonlinear interactions in incompressible magnetohydro- 
dynamics (MHD) are weakened in the field-parallel di- 
rection. The MHD approximation allows to describe 
the large-scale dynamics of astrophysical plasmas, i.e. 
ionized gases, like the interstellar medium or the solar 
corona. Due to the above-mentioned anisotropy, the 
nonlinear energy transfer in MHD turbulence proceeds 
preferably to larger perpendicular spatial wavenumbers 
[l|-|4|]. Direct numerical simulations (DNS) show that 
the field-perpendicular energy spectrum exhibits self- 
similar inertial-range scaling in wavenumber ~ k~ m with 
to — 5/3 @, i] for weak to moderate B , or to — 3/2 
for strong B . 

Iroshnikov [111 and Kraichnan .12] proposed the first 
theory of the effect of a mean magnetic field on incom- 
pressible MHD turbulence. They remarked that any flow 
can be decomposed into a sum of weakly interacting 
waves with different wavevectors k, the term weak inter- 
action meaning that the characteristic time of deforma- 
tion of the waves is much longer than their periods. This 
led to the prediction of a slow cascade, with a spectral 
slope to = 3/2 different from the Kolmogorov prediction 
to = 5/3. 

This theory used an isotropic measure of the propaga- 
tion time based on the modulus of the wave vector, 

< u >^=< (k.Bo)- 1 >~ fcfi- 1 (1) 

ignoring deliberately the waves with wave vectors per- 
pendicular (or almost perpendicular) to the mean field, 
for which the deformation time should clearly be smaller 
than their period, and hence the interaction strong. This 



was criticized by Goldreich and Sridhar [13|, who denied 
the possibility of the previous weak cascade to occur, and 
argued that the perpendicular strong cascade leading to 

— 5/3 

k ± ' should be the only one present. For a large enough 
mean field, the perpendicular cascade should thus be re- 
stricted to a thin subset around the k± axis in Fourier 
space, the subset becoming thinner when the mean field 
increases. 

More precisely, the subdomain in (fcii, k±) space where 
the perpendicular cascade is believed to occur is defined 
by a critical balance [3] between the characteristic time 
of nonlinear interaction, t^l — (k±u\)~ 1 , and the Alfven 
time ta — (&||-0o) _1 , the fluctuations becoming corre- 
lated along the guide field up to a distance ~ B tnl, 
where u\ is the typical magnitude of fluctuations at the 
scale A ~ l/k±. Assuming a scaling law in the perpendic- 
ular direction, and spectral transfer dominated by strong 
coupling, i.e., x = Ta/tnx ^ 1> one obtains the 3D spec- 
trum: 

E 3 (k h k ± ) = k- m -i- 1 f( X ) (2) 

where \ = k ~ q kj_k7^b rm8 / Bq, with to = 5/3 and q = 
2/3, f{x) 1 for |x| > 1 and f( X ) negligible for \ X \ < 1. 

Eq. ([2]) actually suffers from two limitations when the 
mean field Bq is significantly larger than the magnetic 
fluctuation rms value. First, the spectral slope is ob- 
served to become to = 3/2 [7H10( instead of the strong 
cascade value m = 5/3; second, the time-scales ratio 
X = ta/tjsh, becomes significantly smaller than unity 
[14j in the excited part of the spectrum, showing that 
the strict critical balance condition ^ = 1 is too restric- 
tive to describe the anisotropy of the cascade, or, in other 
words, that the cascade is more extended in the oblique 
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directions than predicted by the critical balance condi- 
tion. This has led several authors to suggest modifica- 
tions which either still assume that the anisotropy is dic- 
tated by the critical balance condition [l5|, 16 1, or propose 
that the time-scales ratio x decreases with increasing Bq 
17| . All these phenomenologies predict a spectral form 
different from Eq. ([2]), with different spectral slopes in 
the perpendicular and parallel directions. 

In the present paper, we analyze the angular spectrum 
and find by taking slices along the radial directions that 
a unique spectral form (and slope) holds in all directions, 
with the radial power-law range extent depending on the 
angle with the mean field. As a result, a significant por- 
tion of this spectrum lies in a domain where the time- 
scales ratio x is sub-critical, that is, much smaller than 
unity. 

This study focuses on representative states of fully- 
developed turbulence permeated by a strong mean mag- 
netic field with Bq = 5 6 rms from high-resolution di- 
rect numerical simulations of quasi-stationary MHD tur- 
bulence forced at large scales. The forcing is realized 
by freezing all modes (velocity and magnetic field) with 
k < 2 in an energetically roughly isotropic and equipar- 
titioned state. The driving of magnetic energy could be 
realized physically by large-scale fluctuations of electri- 
cal current although here it is mainly applied to achieve 
a state of approximate equipartion of kinetic and mag- 
netic energy. Decaying test simulations have confirmed 
that the forcing does not modify the results presented in 
the following. The dimensionless equations of resistive 
MHD formulated with the vorticity w = Vxv and the 
magnetic field b are given by 

d t ui = V x [v x lj — b x (V x b)] + pA ui , 
d t b = V x (v x b) + ijAh , 
V-v = V • b = 0. 

The equations are solved by a standard pseudospectral 
method with spherical mode truncation to alleviate alias- 
ing errors. The numerical resolution is 1024 2 x 256 col- 
locations points with reduced resolution in the direction 
of Bo Q and with kinematic viscosity /i and resistiv- 
ity r] set to (i = ?y = 9 • 10~ 5 . The analyzed data 
is the temporal average of five snapshots of the three- 
dimensional Fourier energy distribution taken equidis- 
tantly within about four to five field-perpendicular large- 
scale turnover times, To j_, of quasi-stationary turbulence 
where T , 7 = LoAVms = VK> 3/2 / dk'<5(fc;)K(k')| 2 , 
7 € {x,y,z}, cf. pj], and w rms « b rms = 1 with 
Tq_±_ ~ 1.5, Tq J ~ 1.7. The normalized cross helicity 

p= {■vh)/{(v 2 ) 1 l 2 {b 2 ) 1 ' 2 ) and the Alfven ratio (v 2 )/(b 2 ) 
fluctuate around 23% and 93%. 

The three-dimensional (3D) energy spec- 
trum, E 3 (k x , k y , k z ), relates to the total energy 
E tot = Jd 3 x(v 2 +b 2 )/2 = Jd 3 kE 3 (k x ,k y ,k z ). 

Fig. [T] shows contour levels of £3 in two mutually 





(b) k,= 



100 200 300 400 500 



100 200 300 400 500 
k, 




100 200 300 400 500 
k, 



FIG. 1: Energy contour levels of three-dimensional spectral 
energy density Ez(k x ,k y ,k z ) : (a)plane k z = 0, (b) plane 
k y — 0, (c) average of energy density over all planes containing 
Bo. (d) anisotropic scaling law between wave numbers k± 
and A; 1 1 (see text), with two compensated scalings: k\\(k±)kj_ 
(diamonds) and ku (k±)k^ 2 ^ 3 (crosses). 



orthogonal planes containing the origin. The field- 
perpendicular k x -k y plane (Fig.[TJa) displays an isotropic 
energy distribution, as expected. Anisotropy induced 
by Bo appears in planes containing the Bq direction 
(Fig. Eh). 

Spectral anisotropy is traditionally diagnosed by ID 
spectra, e.g., E±(k x ) = / d 3 k' E 3 5(\k x \ — k' x ) and 
E\\(k z ) = f d 3 k'E 3 S(\k z \ — k' z ). These are shown com- 
pensated by fc 3 / 2 (a) and by fc 5 / 3 (b) in Fig. The per- 
pendicular spectrum exhibits a power-law with m = 3/2, 
in agreement with previous works , while the paral- 

lel spectrum does not show any convincing scaling range. 
An anisotropic scaling law is build from the previous ID 
spectra by plotting (Fig.[TJd) the modes (k±,k\\) sharing 
the same ID energy density [lj|. The anisotropy expo- 
nent q is seen to lie between q = 1 and q = 2/3, which is 
also obtained in [3] for B = 5. 

While planar integration yields some information 
about anisotropy, it mixes all wavenumbers perpendic- 
ular to the chosen direction and thus blurs the separa- 
tion between inertial and dissipative scales if the dissi- 
pative scale is not constant over the planar domain of 
integration. More importantly, no information on in- 
termediate directions between parallel and perpendic- 
ular is available. Thus, spherical coordinates (fc, 9, <fi) 
with respect to the mean field axis along e z are con- 
sidered. As E 3 is isotropic in the azimuthal plane, 
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the 0-dependence of E 3 (k,9,cj)) is eliminated by aver- 
aging over <f> € [0, 27r] which strongly decreases sta- 
tistical noise and yields the (/)-averaged 3D spectrum 
E 3 (k,9) = 1/(271") J dcf>E 3 (k,9, <fi) whose isocontours in 
are shown in Fig. [TJc. We define the cor- 
responding one-dimensional (ID) spectrum E(k, 9) as 
E(k,6) = k 2 E 3 {k,9). The total energy is thus E tot = 
2ir J k 2 dk f* /2 E 3 (k,6) sm(6)d6. 





FIG. 2: Plane-integrated one-dimensional perpendicular and 
parallel energy spectra compensated by (a) k 3 ^ 2 and (b) k'^ 3 , 
respectively. The symbol k stands for the respective field- 
perpendicular and field-parallel wavenumber 



The properties of E 3 (k, 9) are shown in Fig. [3] A scal- 
ing range is seen starting at ko ~ 4-8 down to a dissipa- 
tive wavenumber kd(9). The ID scaling exponent m 



E 3 (k,9) = A(9)k 



-ra-1 



(3) 



is independent of 9. The anisotropy appears at fixed k 
as a ^-dependence of the spectral amplitude A(6) and 
of the dissipative wavenumber kd{9). Normalizing the 
wavenumber by kd shows that the spectrum follows a 
single functional form F whatever 9: 



E 3 (k,0) = F(k/k d ) = F {k/k d {9)) 



-m-2 



(4) 



where Fq is a constant (the amplitude of the spectrum at 
the dissipative scale). Note that the first equality holds 
also beyond the dissipative range, the second being valid 
for k < kd- A corollary is that : 



A(9) oc k d {9) 



m+2 



(5) 



Fig. [3la shows the self-similar wavenumber intervals of 
all E 3 (k,9) spectra, starting in the range ko ~ 4-8 with 
a slight dependence on 9. The dissipative wavenumber 
kd{9) is estimated by locating the maximum of k 2 E{k) in 
each direction 9: varying 9 from 7r/2 to leads to a drop 
of kd{&) from about 100 to 14 while the spectral energy 
at fixed k decreases in the inertial range from 1 to 10~ 3 . 
Fig. [3lb indicates that the spectral ^-dependence can be 
nearly eliminated by normalizing with kd (Eq. ([3])). The 
relation between spectral amplitude A(9) and kd{9) as 



given by Eq. ([5]) is confirmed by Fig.[3lc which shows that 
k 7 J 2 (dotted curve) closely follows A(9) = E 3 (k,9)k 7 / 2 
(Eq. ©), with k < 20 to eliminate the dissipative range 
in the parallel direction. 

A simple model of anisotropic spectrum with a spec- 
tral exponent being the same in all directions has been 
proposed previously in the context of shell models of tur- 
bulence I19H as well as solar wind turbulence 
reads: 



A{9) = {cos 2 9/e 2 + sin 2 9)- {l+m ^ 



(6) 



We tried to use this model to adjust the energy con- 
tours of our numerical simulations, and found that the 
global anisotropy between the perpendicular and parallel 
amplitude requires e = 0.158. However, the model fails 
to reproduce correctly the detailed angular anisotropy, 
that is, the contours in oblique directions, because our 
energy contours differ much from ellipsoids, as is seen in 
Fig. [3jd which shows a zoom of the energy contours as 
solid lines (the model of Eq. ([6]) would produce circular 
contours in this figure). However we found that switch- 
ing from the 2 to 3 for the exponents of the sine and 
cosine as: 



A{6) = {cos 3 9/e 2 + sin z 9)- {l+m W 



(7) 



leads to a reasonable good fit to the simulation results 
in all directions, as seen both in the dotted-dashed curve 
in Fig. [21c for the amplitude variation vs 9 and in the 
dotted contours in the (ku,k±) plane in Fig. Od. 

The energy contours of the critical balance spectrum 
(Eq. [21 with B — 5 and fc = 5) are also represented by 
dashed lines in Fig. Eld. They isolate a small cone about 
the fen axis in the whole plane, so excluding a large part 
of the angular structure of the true angular spectrum. 

Note that the dissipative wavenumber is determined up 
to an error of about a factor two (due to errors in interpo- 
lating the spectrum) , which leads to the noisy appearance 
of the curve of kd in Fig. [31c, to the finite thickness of the 
normalized spectra in Fig. [31b and to variations of about 
a factor 4 in the constant Fo in Eq. (J4|. 

To determine the scaling exponent m of E(k, 9) ~ 
k~ m , the ID spectra averaged over four 9- intervals and 
compensated by fc 3 / 2 and fc 5 / 3 are shown in Fig. [4j The 
spectrum with 9 = tt/2 is represented by the bold line. 
The m = 3/2-scaling is seen to be dominant, except 
possibly for group B which follows m = 5/3 (Fig. [4ja). 
The extent of the 3/2 inertial range is shown by oblique 
dashed lines; it is bounded on the right by kd, and on the 
left by an intermediate range which separates the inertial 
from the forcing range. The start of the inertial range is 
thus growing from k ~ 5 to 10 for 9 — > tt/2. 

Increasing the mean field up to Bq = 5\/2 in test 
simulations (not shown) leads to further decrease of 
the power-law range in the parallel direction, with the 
perpendicular range increasing slightly, and the parallel 
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FIG. 3: Details of spectral properties: (a) Ez(k,9) for 6 
ranging from to tt/2, (b) E 3 (k/k d ,6) (c) k 7 J 2 (dotted), 
k 7/2 E s (k,6) for 8 < k < 20 (solid), Eq.© with m = 3/2 
(dotted-dashed) (d) energy contour levels of E%(kn, k±): sim- 
ulation data (solid); Eq. [5] (dashed line, the oblique line trac- 
ing the boundary x=l with ko — 5, Bo = 5); Eqs. I3I7I with 
P = 3,e= 0.158 (dotted) 





e ~ l/i?o- Increasing again Bo up to 10 confirmed this 
trend, but the parallel range becomes too small in that 
case to allow a good determination of the associated dis- 
sipative wavenumber. This precludes checking that the 
spectrum normalized by k d is angle- independent for small 
8. This difficulty could be alleviated by increasing the 
numerical resolution which would allow increasing the 
Reynolds number. We choose instead below to compare 
with ideal MHD simulations with 512 3 resolution and 
with Bq = 5 and 10. 
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FIG. 5: Anisotropy in ideal runs; (a-b) Bo = 5 (c-d) Bo = 10. 

(a) -(c): fc" /3 (dotted), k 11/3 E 3 (k,9) (solid) (a) for 6 < k < 12 
and (b) for 5 < k < 8, Eq.© (dashed-dotted) with m = 5/3, 
and (a) e = 0.158 (c) e = 0.158/2. (b)-(d): energy contour 
levels of Es(k^,k±): simulation data (solid); Eq. [2] (dashed 
line, the oblique line tracing the boundary x = \ with ko = 5, 
Bo = 5 (b) and B = 10 (d)); Eqs. Ed with /3 = 3, e = 0.158 

(b) and e = 0.158/2 (dotted) 



FIG. 4: Spectra averaged in four subsets of directions: A: 
14° < 8 < 39°; B: 42° < 9 < 65°; C: 67° < 6 < 76°; D: 
79° < 6 < 90°; direction 6 = 90° as thick curve, (a): spectra 
compensated by k 3 ^ 2 , the two oblique dashed lines indicating 
roughly the inertial range; (b): spectra compensated by k 5 ^ 3 



range decreasing substantially, so that the ratio of both 
ranges varies with Bo as 

k d (7r/2)/k d (0) = (A(n/2)/A(0))^ ~ B Q (8) 

Note that the fit by Eq. ([7]) remains as good as in Fig. [3jd 
after decreasing e by a factor \/2, as expected since 
Eq. implies k d (n/2)/k d (0) = l/e hence from Eq. © 



In the ideal case, Fourier space is divided in a large 
scale range presenting spectral properties close to those 
of a standard turbulent spectrum with dissipation, and 
a small scale range where the spectral slope increases, 
the latter scales playing the role of a dissipative range 
(cf. [22I ] in the hydrodynamic case). The boundary be- 
tween the two domains slowly shifts with time to ever 
larger scales. It can be identified with the dissipative 
wavenumber, thus determined here as the minimum of 
the ID spectrum. The simulations are initialized with a 
quasi-stationary state of the resistive run and are con- 
tinued without any dissipation and with the chosen Bq 
in the same numerical setup until the energetically rising 
small scales begin to pollute the scaling region. Choosing 
the appropriate time, power-law ranges in all directions 
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can be properly identified even with a large field Bq = 10. 
The resulting power-law is found to be now m = 5/3 
both with B Q — 5 and f 0, contrary to the resistive runs. 
This difference in scaling between the resistive and ideal 
runs might be attributed to a different role played by the 
bottlenkeck effect [23[ in these two setups. All reported 
findings, cf. Eqs.([3][8]), are however confirmed when set- 
ting m = 5/3, as seen in Fig. [5] which shows (cf. Fig. (He) 
the anisotropy in two ideal runs with Bq = 5 (a) and 
Bo = 10 (b). 

Let us come back to the resistive case. As already 
found in [3], we find that the excited part of the kn,k± 
space is not restricted to regions where x > 1. The 
important point is that the 3D-energy contours, and as 
well the boundary of the power-law range, ignore the iso- 
contours of x as shown in Fig. |3] (see also Fig[5]3,c): the 
form of the 3D-energy contours, as well as the angle- 
independent spectral slope, actually suggest an isotropic 
cascade, that is, a cascade along radial directions. To 
test this idea, we define a ^-dependent effective Reynolds 
Re a oc A(0) 1/2 based on the energy density A(9)kg m ^ 2 
at fco- For 9 = — s- it/ 2, Reo increases by a factor 30, 
while kd grows by a factor fO, as 

k d oc Re^ (9) 

where a — 2/(m + 2) ~ 1/2, m being the ID slope. 
The exponent in Eq. © is substantially smaller than the 
value a — 3/4 for m — 5/3 (or 2/3 in the case m — 
3/2) obtained by equating the input flux at k and the 
dissipative flux e ~ i/k^u(kd) 2 ). This means that if 6 
increases from zero to ir/2, the inertial range increases 
more slowly than it would if the dissipative loss at small 
scales would balance the input energy rate at the (fco — 8) 
wavenumber which marks the large-scale boundary of the 
inertial range. Hence, for 8 —> ir/2 the nonlinear radial 
energy flux must be depleted while the contrary is true 
for the parallel directions. 

By examining the solar wind turbulence, it has been 
shown 20] that the power-law index of the fluctuation 
spectra is independent of the angle between the wave vec- 
tor and the mean interplanetary magnetic field, which 
is fully compatible with the results reported here from 
direct simulations. Other studies of solar wind turbu- 
lence have however reached a different conclusion [2~ij |. 
The latter study used wavelet transforms, which allows 
to define parallel and perpendicular directions with re- 
spect to local averages of the magnetic field. Indeed, it 
has been argued by [H, 0] that the critical balance phe- 
nomenon, and the associated spectral laws, (in particular 
the anisotropy index q = 2/3 relating the perpendicular 
and parallel wave numbers) emerge only when consider- 
ing, instead the mean field, the local average field. In 
the work by [13] who consider as we do only the global 
mean field, the q = 2/3 law appears clearly only when 
the mean field is large enough, which can be explained 



by the fact that in this limit the local and mean field ap- 
proach coincide. According to this viewpoint, the results 
reported here should be a simple artifact of the fact that 
our frame is not attached to the local mean field, but to 
the global mean field. 

However, while the effect could indeed appear for inter- 
planetary turbulence where the fluctuation level is high, 
it is hardly the case here, since Bo/b rms = 5 or 10. Be- 
sides, it is remarkable that the fit by our anisotropy func- 
tion in Fig. [S] is as good when Bq — 5b rms as when 
Bq = 10b rms which is an indication that the spectral 
properties of the turbulence are correctly revealed by us- 
ing our method. 

A possible way to reconcile both pictures is the fol- 
lowing. It takes into account the fact that there is a 
second parameter, the Reynolds number. Indeed, as 
the anisotropy increases with Bq, we find that the scal- 
ing range in the parallel direction decreases accordingly 
(Eq. (j8])), which implies that, at a fixed viscosity, the par- 
allel power-law range disappears at even moderate Bq. 
In our case, kd± — 100 and fc^n — 10 when Bq — 5b rrns , 
while the start of the scaling range is fco — 5 — 8, pre- 
venting to increase significantly Bq. The regime at high 
Bq will thus depend on the Reynolds number [Re). In- 
creasing Bq at fixed perpendicular Re depletes the field- 
parallel cascade and could ultimately lead to critically 
balanced turbulence. If -Bo Re are however large 
enough, e.g. under astrophysical conditions, allowing for 
scaling in all directions then the properties described 
in this work are expected to hold. We thus propose 
direction-independent scaling for high Re and critically 
balanced turbulence at low Re. The crossover Re- value is 
expected to increase with Bq . A strong indication in this 
sense is found in [3] where the anisotropy scaling expo- 
nent q is seen to be between 1 and 2/3 at Bq = 5b rms as 
here (Fig. [3d), while it begins to cluster around 2/3 at 
BQ/b rms > 10. Such a good agreement with q = 2/3 is 
found already at Bo/b rms ~ 1 in [5j because of a moder- 
ate Reynolds number. 

We have reported here two previously unknown prop- 
erties of the MHD angular energy spectra: (i) its func- 
tional form is A(9)f(k) (ii) the anisotropy function A{9) 
is best expressed (Eq. [8|) as the ratio of the perpendicu- 
lar over the parallel power-law range extent, which scales 
linearly with Bq in the Bo/b rms = 5, 10 interval consid- 
ered here. These results offer important tests for future 
theories of anisotropic turbulence. 

We thank G. Belmont, J. Leorat, and A. Busse for 
several fruitful discussions. 
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